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Abstract. Measurement of jet fragmentation property in p + p collisions provides a crucial 
O ' baseline for the study of a possible, modified fragmentation behavior at the presence of the 

2 ■ quark-gluon plasma. We describe the first measurement of jet fragmentation functions in 

■v/s = 200 GeV p + p collisions at RHIC. This is facilitated by extending the linear least square 
unfolding technique employed by high energy and nuclear experiments to multidimensional 
spectra, which allows us to correct the inefficiencies in the jet energy response encountered at 
detectors such as PHENIX. 
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•/^ ' 1. Introduction 

Hadronization of jet production into final state particles, together with the production of 
hard scattering, are complementary parts of the QCD description of particle production. 
Fragmentation function in p + p collisions further provides a baseline for any measurement of 
the fragmentation property of reconstructed jets in heavy ion collisions. The narrow detector 
^ I aperture and the need to operate within the large multiplicity heavy ion background has 

j^ ■ traditionally limited the application of direct jet reconstruction at PHENIX. Fragmentation 

properties of jet production therefore been only measured via two particle correlations by 
PHENIX [H [2]. The Gaussian filter jet reconstruction algorithm [3] provides a cone-like 
algorithm that at the same time reduces the sensitivity to large angle fragments (or the lack 
of such, due to acceptance limit), and is therefore well-suited for both p + p and heavy ion 
measurement at PHENIX. 

Measurement of fragmentation function using a detector with inefficiency due to the lack 
of hadronic calorimetry is difficult, since the full jet energy is not directly accessible, but has 
to be reconstructed from combining track momenta with the electromagnetic cluster energies. 
While the energy carried by long-lived neutral hadrons such as n, K^ are fully lost, the inherent 
momentum resolution and background in a typical tracking system limits the ability to accurately 
measure tracks with p^ 20 GeV /c. The resulting large difference in the true and measured jet 
energy scale is difhcult to correct multiplicatively, but can be addressed by using unfolding 
the measured spectra with the true-to-detector-level jet energy transfer function. This is 
particularly difficult for spectra that contain additional dimensions than the jet energy, such 
as the fragmentation function. As such, none of the RHIC experiment, which all lacks hadronic 
calorimetry, has attempted the measurement of fragmentation functions so far. 



Jet reconstruction using the crfiiter = 0.3 Gaussian filter has been applied to PHENIX for 
the measurement of the inclusive jet cross section in p + p and jet yield in Cu + Cu [H [5]. In 
this proceeding, we report the development of a multidimensional unfolding technique for the 
energy scale correction of fragmentation functions, and the measurement oi p + p fragmentation 
functions for charged particles and neutral electromagnetic clusters using the PHENIX detector. 

2. Jet reconstruction by Gaussian filtering 

The iterative cone algorithm with the size parameter R is known to be equivalent to a local 
maximization of a filter output in (r/, (p) with a special choice of the angular weight function 
to be the flat k^r"^) = 6{R^ — r^) and r'^ = t]'^ + cfp' (note that unlike k(r'^), the filter kernel is 
/i(r^) oc — j dr'^k{r'^) oc max(0, 1 — r'^/R^), and not flat) [6l[7]. Similarly, the Gaussian filter uses 
a Gaussian distributed weighting. The Gaussian weighting takes advantage the property of high- 
Pt jets being collimated emissions of particles, and enhances the center signal to the periphery, 
which is more likely to be contaminated by the event background. By avoiding an sharp radial 
cut-off, the algorithm also becomes analytically collinear and infrared safe (we further verified 
the practical infrared safety using a procedure analogously to [8]). 

\n.p-\- p events and without a large event background, the event transverse momentum density 
is the sum of point-like final state particles pT,i 

PT{v,(t>) = '^PT,iS{v-Vi)S{4>-4>i)- (1) 

The Gaussian filtering of py is the linear-circular convolution of pT{r],<p) with a Gaussian 
distribution 

p^Xr,, <A) = // d^'^'pHV, 'A') exp ( J^-^')l + ^'i^-'t^'n . (2) 
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The output of the filter for a given (77, 0) position is the Gaussian weighted transverse momentum 
in the said event. The local maxima in p|J*(7/, (/>) are the reconstructed jets using the Gaussian 
filter. More detailed discussion of the property of Gaussian filter for jet reconstruction can be 
found in [21 [5] . 

3. Multidimensional unfolding 

The true particle level and the measured energy scale in the PHENIX central arm differ due to 
inefficiencies, the lack of hadronic calorimetry, and the loss of out-of-acceptance fragments due 
to the A77 = 0.7 pseudorapidity coverage and also the Ac/; = vr partial azimuthal coverage. This 
difference means that a full unfolding is needed to extract the true energy instead of a first order 
multiplicative correction. 

Linear least square unfolding with the Phillips-Tikhonov regularization [9l [10] has been 
widely used by e.g. GURU [llj for the unfolding of spectra. Binned spectra can be mathematically 
regarded as a vector of individual bins, and the linear bin-to-bin migration due to differing energy 
scales is then conveniently expressed as a matrix. For a typical transfer matrix A generated 
using Monte Carlo detector simulation, if the measured spectrum b with the covariance matrix 
C is unfolded using the simple linear least square relation 



minimize || Ax — b||Q, (3) 



(||u||c = u'^C~^u denotes the Mahalanobis distance) the propagation of fluctuations in A would 
result in large, non-statistical fluctuation in the unfolding result x. The regularized unfolding 
determines the unfolded spectrum x by solving the minimization problem 



minimize ||Ax — bH^ -|- r||Lx|| , (4) 



where r is the regularization parameter and L is a hnear operator that describes the amount of 
discontinuity or "noisiness" of the unfolding result. 

Spectra in particle and nuclear physics can cross multiple orders of magnitude and a 
homogeneous L would cause the large magnitude part of the spectra to dominate the 
singular vectors, thus exhausting most degrees of freedom purely to reproduce the variation 
in magnitudes. This can be solved by either scaling L to approximately match the variations in 
X, or to prescale x. We follow the approach used by GURU to prescale the unknowns by 

Xi^ Xi = -j^, (5) 

and the linear system A is therefore inversely scaled by 

fljj I— 7- ttij = CLijX^ [b) 

(column- wise scaling), with the purpose of preventing exceedingly large or small Xi to appear 
numerically. 

The typical choice for L is a second order finite-difference matrix 
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The minimization with respect to the second order derivative describes a continuity constraint 
that restricts the shape of x to be cubic spline-like. The simplest Z)-dimensional generalization 
of L*-^' are the isotropic axial derivatives d'^/dxf.,k = 1,...,Z), and the minimization with 
respect to it behaves similarly to the D-dimensional cubic tensor splines. We found this type of 
regularization to be sufficient for the purpose of unfolding multidimensional spectra with Poisson 
statistics, other choices used e.g. in image restoration are reviewed in [12]. 

The solution to ^ can be written as a least square problem with the matrix pencil (A-|--y/rL), 
and can be solved by a SVD of (A"1A)L-i = USV^, where C = AA^ is the Cholesky 
decomposition of the measurement covariance, and provided L is nonsingular. The software 
package GURU provides an implementation of this method for diagonal C and using L^^^ (which 
is strictly speaking singular, but the rectangular nature allows it to be approximately inverted 
by perturbing L^^' i— t- L'^-* -|- el, |e| <^ 1). However, a D-dimensional continuity constraint for 
A^ measurements consists of at least DN axial derivatives and therefore L becomes a DN x N 
matrix, it is therefore generally not possible to invert L and solve the multidimensional unfolding 
in the fashion of GURU. 

The least square problem with a matrix pencil (A~ A -|- \'^'L) can also be solved by 
the generalized SVD (GSVD) [13] of a matrix pair GSVD(A^ A,L), which simultaneously 
decomposes both matrices into 

A~iA = U(0 SOX"^ 

(8) 

L = v(o s;2)x-^ 

GSVD(A^ A, L) is therefore closely related to SVD((A^ A)L^^), which can be immediately 
seen by comparing the GSVD of the quotient matrix (A~ A)L^^ = U(5]i5]2^ )V^ (L^ being the 
(A~"'^A)-weighted pseudoinverse of L [14j) to the ordinary SVD form (A^"'^A)L^-'^ = USV-^. 
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Figure 1. The PHENIX central arm 
detectors for RHIC Run-5 (year 2004/2005), 
viewed along the beam axis from the south 
towards north. Dark regions indicate 

detectors used for the jet reconstruction: The 
drift chamber (DC), the pad chamber layers 1 
and 3 (PC1/PC3), the ring-imaging Cerenkov 
detector (RICH), and the Pb scintillator 
(PbSc) and Pb glass (PbGl) electromagnetic 
calorimeters. 



Figure 2. The PHENIX jet P{p^\pt) 
transfer matrix for ^/s = 200 GeV and a = 
0.3 Gaussian filter, derived from the geant 
simulation of ss 1.6 x 10^ pythia events. The 
p^ < pt region is dominated by n, K^ energy 
loss. 



With the generalized singular value pairs Si = diag(aj), S2 = diag(/3j), the SVD singular 
values of (A~ A)LT are 7^ = ai//3i. The solution then can be written as the Tikhonov filter for 
the singular values 
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and the mean and covariance matrix of the unfolding solution is 



X = X5]'-'U^(A-Ib) 
Covx = XS'-\s'-^)^X^. 



(9) 



(10) 



We implemented this multidimensional generalization to guru's regularized linear least 
square unfolding based on lapack's dggsvd [15J, which implements the dense GSVD algorithm 
by Bai, Demmel & Zha [Ml [H] . 



4. Experimental setup 

Figure [1] shows the Run-5 PHENIX "central arm" configuration for RHIC Run-5 (year 
2004/2005). PHENIX has two mid-rapidity spectrometers with an aperture of |r/| < 0.35 and 
A(/) = 7r/2 each. The components used for jet reconstruction are the drift chamber (DC), the 
pixel pad chamber layers 1 and 3 (PC1/PC3), the ring-imaging Cerenkov detector (RICH), and 
the electromagnetic calorimeters (EMCal). For the data presented in this paper, DC/PC1/PC3 
provide momentum measurement for charged particles, and the EMCal the energy for photons. 
Pattern recognition and momentum reconstruction in the tracking system formed by DC 
and PCl/3 is performed using the combinatorial Hough transform, with the momentum scale 
determined by the time-of-flight measurement of identified vr^, if^, and pjp- The momentum 



resolution of the tracking system is determined as 5p/p = 0.7% © 1.0%p/(GeV/c). Two 
calorimeter technologies were used, six of the total eight sectors are covered by Pb-scintillator 
(PbSc) calorimeters with resolution of (Te/E = 8.1%/\rE©2.1% and a granularity of Ar/ x Ac/) ~ 
0.01 X 0.01, and two sectors by Pb-glass (PbGl) calorimeters with (Te/E = 5.9%/\/!E©0.8% and 
Ary X Ai;^ ~ 0.008 x 0.008. The intrinsic timing resolution for 1 GeV vr are about 200-300 ps 
for both technologies. 

Since PHENIX currently lacks an in-field tracking capability, conversion electrons in the DC 
can produce a displaced track that has the appearance a high pT track originating from the 
event vertex. We use the information from the RICH and dE/dx measurement to identify and 
remove these tracks. To provide additional suppression at the cross section level of jets with 
Pt > 20 GeV, we use the fact that conversion electrons are geometrically unlikely to coincide 
with the direction of the jet production, and require the reconstructed jet to have a minimum 
multiplicity of 3 particles measured within a radial angle of 60° , and the charged fraction of the 
jet Pt to be below 0.9 to remove single track events, or when the jet is dominated by a large pT 
track. 

The absolute energy scale of the calorimeter clusters are determined using both the 
reconstructed vr*^ masses from the observed vr*^ — )• 77 decays, and checked using the E/p from 
energies from cluster matching RICH identified e^ tracks. Shower shape cuts are used to remove 
clusters. The residual uncertainty in the energy scale is ±3% (syst.). Since the measurement 
extends to very low cross sections, processes such as upward beam interaction can deposit energy 
in the EMCal that randomly coincides with an event. This is suppressed by measuring the time- 
of-flight of the clusters and rejecting those that are out of synchronization with the collision. 

The PHENIX minimum bias (MB) trigger is defined by the coincident firing of the two beam- 
beam counters (BBC) located at 3.0 < rj < 3.9. The Van de Meer/vernier scan method is used 
to measure the BBC cross section, with ctbbc = 22.9 it 2.3 mb (syst.). The efficiency of BBC 
MB trigger on an event containing a jet with p^r^^ > 2 GeV/c is ebbc = 0.86 it 0.05 (syst.) and 
constant with respect to p^^ within that uncertainty. We require the collision vertex to be 
within \z\ < 25 cm along the beam axis, derived from the timing difference between the firing of 
the two BBC. 

5. Inclusive jet cross section in p -\- p 

The data presented in the following sections were obtained from the PHENIX p + p dataset 
from the RHIC Run-5 (year 2004/2005). After removal of bad quahty runs, a total of 1.47 x 10^ 
minimum bias p + p and 1.16 x 10^ triggered p + p events are being used. 

The measurement of inclusive jet cross section in Run-5 p + p collisions uses the combined 
minimum bias and triggered data set, and is based on the regularized unfolding of the 
jet spectrum dNjet/dpl^ '^'^^ by the jet energy scale transfer function P(p^y '^'^'^\p!^ ), which is 
evaluated by geant simulation using pythia 6.4.20 events. A total of 1.6 x lO'^ events 
were simulated with 14 different minimum Q^ settings varying betweeny^Q^ > o.5 GeV /c and 
\/Q^ > 64GeV/c. The resulting transfer matrix Pip^lpT) is shown in figure [21 

The unfolding is equivalent to the inversion of the ID Fredholm equation 

^^-^ ''dj^^'Pij^^'ni^^')^^. (11) 



7 iet,rec / '^ ^'^ '-''jl 

dp^rp J apTrp 

Then the invariant cross section can be evaluated as 

Ed^ai""^ 1 dVJ^* CTBBC 1 1 dNi"^ 



(12) 



dp3 2TTPT dp^^dy A eBBC pj^* iVevt dj^^^ 

where dN/dpT is the unfolding result from (|11|) and 

^ = 2(A7/-2(i)(A(/>/2-2d) (13) 
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Figure 3. PHENIX Run-5 p + p at 

y^ = 200 GeV invariant jet cross section 
spectrum as a function of pT- The 

shaded box to the left indicates the overall 
normalization systematic uncertainty, shaded 
boxes associated with data points indicate 
point-to-point systematic uncertainties, and 
error bars indicate statistical uncertainties. 



Figure 4. PHENIX Run-5 p + p at ^/s = 
200 GeV invariant jet cross section spectrum 
as a function of pT, with comparison to 
|18j . next-to-leading order calculation from 
[19] . and PYTHIA assuming K = 2.5. The 
shaded box to the left indicates the overall 
normalization systematic uncertainty, shaded 
boxes associated with data points indicate 
point-to-point systematic uncertainties, and 
error bars indicate statistical uncertainties. 



the fiducially reduced PHENIX central arm acceptance area. This measurement is detailed in 

M- 

The choice of regularization parameter r can translate into an uncertainty on the low 
frequency, global shape of a spectrum. We address this by evaluating the systematic uncertainty 
from varying r over the entire meaningful range between ~ 4 degrees of freedom up to the 
Nyquist frequency, D being the unfolding dimension. 

Figure [3] shows the PHENIX preliminary p + p jet spectrum measured using the Gaussian 
filter, plotted in invariant cross sections. The shaded box to the left indicates the overall 
normalization systematic uncertainty, shaded boxes associated with data points indicate point- 
to-point systematic uncertainties, and error bars indicate statistical uncertainties. We show the 
unfolded spectrum out to the pT bin where the nominal yield for the number of sampled events 
reaches the level of 1 jet, namely 60 GeV /c. 

Figure [4] shows the same spectrum as in Figure [3l compared against the spectrum from |18], 
the next-to-leading order (NLO) calculation using the small cone approximation (SCA) [19], 
and the leading order pythia spectrum assuming K = 2.5. The comparison to [ISj and NLO 
SCA involve different jet definitions, a residual difference should be expected, even though for 
Pt > 15 GeV/c it appears to be small between filter and cone jets for the Gaussian size a = 0.3 
used in this analysis. Our spectrum is close to [18] within its prp reach. The spectrum also follows 
approximately the shape of the NLO SCA calculation, and the leading order pythia spectrum, 
a K = 2.5 is assumed. However, a more appropriate comparison would involve Gaussian filter 
based NLO calculations, which we plan to perform in the future. 



6. Fragmentation functions in p-|-p 

We use the Run-5 p + p minimum bias data set for the measurement of the fragmentation 
function. For reconstructed jets, we measure the longitudinal contribution of the fragment 
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Figure 5. PHENIX Run-5 p + p at ^/s = 
200 GeV charged (with electron rejection) 
jet fragmentation function with respect to 
z = ti'mtcrPji/p^'^* and vertically scaled 
by c{ii^^) = 10^, A; = 0,...,3. The 
shaded box to the left indicates the overall 
normalization systematic uncertainty, shaded 
boxes associated with data points indicate 
point-to-point systematic uncertainties, and 
error bars indicate statistical uncertainties. 
Biases from the jet level cuts are fully 
quantified in the systematic uncertainties. 
The blue curve indicates the D{z) = Nz'^{l — 
zY (l + ^) fit to data, while the red curve 
shows the same fit from pythia at pj! = 
15 GeV/c. 



Figure 6. PHENIX Run-5 p + p at 

^/s = 200 GeV (electromagnetic) neutral 
jet fragmentation function with respect to 
z = w^\tcvP\\/p^'^^ and vertically scaled 
by c{j}^^) = 10^=, A; = 0,...,3. The 
shaded box to the left indicates the overall 
normalization systematic uncertainty, shaded 
boxes associated with data points indicate 
point-to-point systematic uncertainties, and 
error bars indicate statistical uncertainties. 
Biases from the jet level cuts are fully 
quantified in the systematic uncertainties. 
The blue curve indicates the D{z) = Nz"'{l — 
z)^ (1+2) fit to data, while the red curve 
shows the same fit from pythia at pj. = 
15 GeV/c. 



momentum to the jet as 
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(14) 



This choice preserves the identity ^ z = 1 for a weighted jet definition, compared to 
traditional algorithms that would assign zero or unitary weights for nonconstituent and 
constituent fragments. For practical purpose, both pu and p!^ are binned logarithmically in 

the dN/{dp\\dpl^ ) distribution, which allows the p\\/pT division to be performed as bin shifts. 
For the narrow mid-pseudorapidity acceptance of the PHENIX central arms, the difference in 
z = pii/p''^* and z = p\\/pj^ division is < 6% and negligible within our uncertainties. 

We use RICH to reject electrons for the measurement of charged fragmentation function, 
since in PHENIX, they are produced predominantly from 7 conversions in the beam pipe and 



Dalitz decays, therefore not associated with the charged fragmentation of jets. At the moment, 
we also do not reconstruct 7 pairs from n^ decays, although we envision to provide identified 
neutral fragmentation function measurements in the future. 

Since the PHENIX detector resolution for single particle far exceeds the resolution for jet 
energy, we assume a diagonal transfer function for particle energies in the transfer function 

pii^T''''',i^nPT,p\\) = PiP^'ni^T^pwmr - p\0^ (i^) 

although the use of a full transfer function is planned for the future. The same GEANT simulation 
result as for the jet spectrum is used to evaluated the so reduced conditional transfer function 

7^/ iet,reci iet \ 

P{Pt \Pt 'P\\)- 

The {pj, ,p||) distribution is then unfolded by the regularized inversion of the 2D equation 

"'-'■^ II ■ iet 7 J-,/ iet.rec reel iet n "'-'■^ 



Similarly, the jet spectrum is unfolded by inverting (jlip . The per-jet normalized fragmentation 
function is then extracted by dividing the unfolded {pj, ,p\\) distribution by the unfolded jet 
spectrum: 



dTYjetT^ clN _ 1 /diVjct\ ^jet dN 



^W = -^l^M^fi^ -if- = -^T^M^fl^ P^ 



e{z,p!^ ) y dp!^ J dpl^ dz e{z,p!^ ) y dp!^ J dp!^ dp\ 



(17) 
=Pii/pj<=* 



where €{z,pj, ) is the single particle efficiency evaluated by geant simulation. While the 
efficiency for neutral clusters are constant at e = 0.81 it 0.04 (syst.), the efficiency for charged 
tracks exhibits a strong z-dependence both due to the magnetic field causing low px tracks to 
be bend out of the azimuthal acceptance, and tracking cuts resulting in the rejection of high p^ 
TT^ misidentified as electrons. For pj, = 18 GeV/c, the efficiency for charged tracks saturates 
for z > 0.3 to e = 0.58 ± 0.06 (syst.). 

Figure [5] shows the p + p at -y/i = 200 GeV charged (with electron rejection) jet fragmentation 
function, vertically scaled by c{p!^ ) = lO'^, A; = 0, . . . , 3. The shaded box to the left indicates the 
overall normalization systematic uncertainty, shaded boxes associated with data points indicate 
point-to-point systematic uncertainties, and error bars indicate statistical uncertainties. Biases 
from the jet level cuts are fully quantified in the systematic uncertainties. The blue curve 
indicates the D{z) = A^z"(l — z)^ (l + ^) fit to data, while the red curve shows the same fit 
from PYTHIA at pi^^ = 15 GeV/c. 

Figure [6] shows the same plot, for the p + p at ^/s = 200 GeV (electromagnetic) neutral jet 
fragmentation function. In both set of fragmentation measurements, the fragmentation function 
for p!^ > 15 GeV /c agree well within our current uncertainty with pythia. We could reach a 
maximum z ~ 0.81 with our measurement. 

7. Summary 

We presented the first measurement of jet fragmentation function in p + p collisions performed 
at RHIC. It demonstrates for the ability of PHENIX to access high-z jet fragmentation property 
for both charged particles and photons. Techniques to unfold multidimensional spectra such as 
the fragmentation function has been discussed, which significantly facilitated us to make such a 
measurement. 

Restricting the data to the minimum bias set currently constrains our statistics and ability 
to lower the unfolding systematic uncertainty. Our current development of the trigger efficiency 
correction for fragmentation function measurement would allow us to address this problem in 
the future. 
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